set.seed(123)
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
knitr::opts_chunk$set(fig.width = 10, fig.height = 10, fig.dpi = 600)
require(EnhancedVolcano)
require(DESeq2)
require(limma)
require(IHW)
require(tidyverse)
require(dplyr)
require(ggplot2)
require(Hmisc)
require(org.Hs.eg.db)
require(clusterProfiler)
SIG_THR <- 0.05
LFC_THR <- 1
MIN_COUNT <- 10
MIN_SUMMED_ROW_COUNT <- 3
# Source: https://rdrr.io/github/ttecon/ttr/src/R/geom_split_violin.R
#
#' Split Violin plot
#'
#' This function is to plot the split violin graph.
#' It is built on [ggplot2::geom_violin()],
#' and adopted from jan-glx's stackoverflow answer \url{https://stackoverflow.com/a/45614547}
#' A violin plot is a compact display of a continuous distribution.
#' It is a blend of [ggplot2::geom_boxplot()] and [ggplot2::geom_density()]:
#' A violin plot is a mirrored density plot displayed in the same way as a boxplot.
#' A split violin plot allows for odd violins to be plotted to the left,
#' and even violins to be plotted to the right
#'
#' @param mapping Set of aesthetic mappings created by `aes()` or `aes_()`.
#' If specified and `inherit.aes = TRUE` (the default),
#' it is combined with the default mapping at the top level of the plot.
#' You must supply mapping if there is no plot mapping.
#' @param data The data to be displayed in this layer. There are three options:
#'
#' If `NULL`, the default, the data is inherited from the plot data
#' as specified in the call to [ggplot2::ggplot()].
#'
#' A `data.frame`, or other object, will override the plot data.
#' All objects will be fortified to produce a data frame.
#' See [ggplot2::fortify()] for which variables will be created.
#'
#' A function will be called with a single argument, the plot data.
#' The return value must be a `data.frame`, and will be used as the layer data.
#' A function can be created from a formula (e.g. `~ head(.x, 10)`).
#' @param stat Use to override the default connection between `geom_density` and `stat_density`.
#' @param position Position adjustment, either as a string,
#' or the result of a call to a position adjustment function.
#' @param ... Other arguments passed on to `layer()`.
#' These are often aesthetics, used to set an aesthetic to a fixed value,
#' like `colour = "red"` or `size = 3`.
#' They may also be parameters to the paired geom/stat.
#' @param draw_quantiles If `not(NULL)` (default), draw horizontal lines
#' at the given quantiles of the density estimate.
#' @param trim If `TRUE `(default), trim the tails of the violins to the range of the data.
#' If `FALSE`, don't trim the tails.
#' @param scale If `"area"` (default), all violins have the same area (before trimming the tails).
#' If `"count"`, areas are scaled proportionally to the number of observations.
#' If `"width"`, all violins have the same maximum width.
#'
#' @param na.rm If `FALSE`, the default, missing values are removed with a warning.
#' If `TRUE`, missing values are silently removed.
#' @param show.legend logical. Should this layer be included in the legends?
#' `NA`, the default, includes if any aesthetics are mapped.
#' `FALSE` never includes, and `TRUE` always includes.
#' It can also be a named logical vector to finely select the aesthetics to display.
#' @param inherit.aes If `FALSE`, overrides the default aesthetics,
#' rather than combining with them.
#' This is most useful for helper functions that define both data and aesthetics
#' and shouldn't inherit behaviour from the default plot specification, e.g. [ggplot2::borders()].
#'
#' @export
geom_split_violin <- function(mapping = NULL,
data = NULL,
stat = "ydensity",
position = "identity",
...,
draw_quantiles = NULL,
trim = TRUE,
scale = "area",
na.rm = FALSE,
show.legend = NA,
inherit.aes = TRUE) {
ggplot2::layer(
data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(
trim = trim, scale = scale, draw_quantiles = draw_quantiles,
na.rm = na.rm, ...
)
)
}
#' @format NULL
#' @usage NULL
GeomSplitViolin <- ggplot2:::ggproto("GeomSplitViolin",
ggplot2::GeomViolin,
draw_group = function(self, data, ..., draw_quantiles = NULL) {
data <- transform(data, xminv = x - violinwidth * (x - xmin),
xmaxv = x + violinwidth * (xmax - x))
grp <- data[1, "group"]
newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv),
if (grp %% 2 == 1) y else -y)
newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
1))
quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")),
drop = FALSE]
aesthetics$alpha <- rep(1, nrow(quantiles))
both <- cbind(quantiles, aesthetics)
quantile_grob <- GeomPath$draw_panel(both, ...)
ggplot2:::ggname("geom_split_violin",
grid::grobTree(GeomPolygon$draw_panel(newdata, ...),
quantile_grob))
}
else {
ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
}
}
)
Feature counts were obtained using Rsubrea::featureCounts(). Strandness was inferred using salmon utilized in istrand script. All samples were inward inversely stranded libraries.
# load count matrix
counts <- read.csv("workdir/feature_counts/feature_counts.csv",
row.names = 1)
cols <- colnames(counts)
cols <- gsub("_sorted.bam", "", cols)
cols <- gsub("X", "", cols)
colnames(counts) <- cols
# change entrezIDs to gene symbols if available
id.to.symb <- bitr(rownames(counts),
fromType = "ENTREZID",
toType = c("SYMBOL"),
OrgDb = org.Hs.eg.db,
drop = FALSE)
gene.symbols <- ifelse(is.na(id.to.symb$SYMBOL),
id.to.symb$ENTREZID,
id.to.symb$SYMBOL)
rownames(counts) <- gene.symbols
counts[1:5, 1:5]
df <- read.csv("workdir/clinical.csv")
df$OS_median <- ifelse(df$OS > median(df$OS), "above", "below") %>%
as.factor() %>% relevel(ref = "above")
df$DFS_median <- ifelse(df$DFS > median(df$DFS), "above", "below") %>%
as.factor() %>% relevel(ref = "above")
df$CR <- df$CR %>% as.factor
str(df)
## 'data.frame': 165 obs. of 17 variables:
## $ ID : chr "8" "11" "30" "37" ...
## $ therapy : chr "PC" "PC" "PC" "PC" ...
## $ age : int 55 65 34 48 46 66 59 64 45 53 ...
## $ CR : Factor w/ 2 levels "no","yes": 1 2 2 2 2 1 1 1 2 2 ...
## $ PS : chr "no" "yes" "yes" "yes" ...
## $ OS : int 104 1193 467 2320 3750 385 790 1138 1847 2351 ...
## $ DFS : int 0 393 226 815 1468 0 0 0 1429 528 ...
## $ relapse : chr NA "yes" "yes" "yes" ...
## $ death : chr "yes" "yes" "yes" "yes" ...
## $ TP53_mm : chr "yes" "yes" "no" "no" ...
## $ p53_acc : chr "yes" "yes" "no" "no" ...
## $ type : chr "HGSOC" "HGSOC" "LGSOC" "LGSOC" ...
## $ microinvasion: chr NA NA NA NA ...
## $ implant : chr NA NA NA NA ...
## $ primary_tumor: chr NA NA NA NA ...
## $ OS_median : Factor w/ 2 levels "above","below": 2 2 2 1 1 2 2 2 1 1 ...
## $ DFS_median : Factor w/ 2 levels "above","below": 2 1 2 1 1 2 2 2 1 1 ...
Because of the difference in expression profiles between HGSOC and other types and the fact that HGSOC is the largest sample group this group was selected for subsequent analysis. Number of HGSOC samples was 140.
dff <- df %>% filter(type == "HGSOC")
nrow(dff)
## [1] 140
Firstly, patients with OS above median were compared to those with OS below median. Low count rows (genes) were filtered out to avoid noise. Exact experiment design is shown below.
dds <- DESeqDataSetFromMatrix(countData = counts[dff$ID],
colData = dff,
design = ~OS_median)
keep <- rowSums(counts(dds) >= MIN_COUNT) >= MIN_SUMMED_ROW_COUNT
dds <- dds[keep, ]
dge <- DESeq(dds)
res <- results(dge, filterFun = ihw)
head(res, 0)
## log2 fold change (MLE): OS median below vs above
## Wald test p-value: OS median below vs above
## DataFrame with 0 rows and 7 columns
Number of significantly differentiating genes and volcano plot are show below.
signif <- res %>% as.data.frame %>% filter(abs(log2FoldChange) > LFC_THR,
padj < SIG_THR)
write.table(signif, file = "workdir/DGEA/prog.tsv", quote = FALSE)
length(signif %>% row.names)
## [1] 219
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
xlab = "LFC",
ylab = "FDR",
pCutoff = SIG_THR,
FCcutoff = LFC_THR,
legendLabels = c("NS", "LFC", "FDR", "FDR and LFC"))
resLFC <- lfcShrink(dge, coef = 2, type = "apeglm")
DESeq2::plotMA(resLFC, ylim = c(-5, 5))
Top 10 significant genes were selected and the differences between 2 groups are plotted below according to (a) log count (b) log SF normalized count (c) VST normalized count.
top_df <- signif %>% arrange(padj) %>% head(10)
top_df
dds <- estimateSizeFactors(dds)
norm.counts <- counts(dds, normalized = TRUE)
vst.counts <- dds %>% vst %>% assay
top_genes <- top_df %>% rownames
top_genes_data <- data.frame()
for (gene in top_genes) {
d <- plotCounts(dds, gene=gene, intgroup="OS_median", returnData=TRUE)
d$ID <- rownames(d)
d$gene <- gene
new_df <- data.frame(gene=d$gene, ID=d$ID,
OS_median=d$OS_median,
count=d$count)
top_genes_data <- rbind(top_genes_data, new_df)
}
top_genes_data$norm_count <- (norm.counts[unique(top_genes_data$gene), ]
%>% t)[dff$ID,] %>% as.vector
top_genes_data$vst_count <- (vst.counts[unique(top_genes_data$gene), ]
%>% t)[dff$ID,] %>% as.vector
ggplot(top_genes_data, aes(x=gene, y=log(count), fill = OS_median)) +
geom_split_violin(trim = 0) +
geom_boxplot(alpha = 0.3, width = 0.2, outliers = FALSE) +
theme_classic() +
ylab("log count") +
guides(fill = guide_legend(title = "OS median")) +
theme(axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank())
ggplot(top_genes_data, aes(x=gene, y=log(norm_count), fill = OS_median)) +
geom_split_violin(trim = 0) +
geom_boxplot(alpha = 0.3, width = 0.2, outliers = FALSE) +
theme_classic() +
ylab("log SF normalized count") +
guides(fill = guide_legend(title = "OS median")) +
theme(axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank())
ggplot(top_genes_data, aes(x=gene, y=vst_count, fill = OS_median)) +
geom_split_violin(trim = 0) +
geom_boxplot(alpha = 0.3, width = 0.2, outliers = FALSE) +
theme_classic() +
guides(fill = guide_legend(title = "OS median")) +
ylab("VST normalized count") +
theme(axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank())